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I. INTRODUCTION 

In two seminal papers Svetitsky and Yaffe have tentatively 
linked the finite-temperature phase transitions in 'hot' gauge 
theories to the simpler order-disorder phase transitions of 
spin models In the general case their conjecture may 

be stated as follows: The effective theory describing finite- 
temperature Yang-Mills theory with gauge group SU{Nc) 
in d + 1 dimensions is a spin model in d dimensions with 
a global symmetry group given by the centre Z{Nc) of the 
gauge group. A somewhat stronger version of the conjecture 
can be formulated if the phase transitions in questions are of 
second order. In this case Yang-Mills theory and spin model 
fall into the same universality class and the critical exponents 
coincide. This has been convincingly demonstrated for SU (2) 
011 md = 3. 

However, continuous phase transitions in hot gauge theories 
are not generic and hence the universality statement is almost 
empty - at least in 3+1 dimensions ^5^. On the other hand, 
the more general version of the statement has already been 
used by Svetitsky and Yaffe to argue that the phase transition 
for SU (3) Yang-Mills theory must be first order for d = 3 as 
there is no Z(3) RG fixed point in this case. Since then this 
has been firmly established by a number of lattice calculations 
008,9, 10, 11]. 

The conjecture implies that the effective theories may be 
formulated as 'Polyakov loop models ' O El El. 

For 

SU (Nc) this means that the macroscopic dynamical variables 
have to reflect the complete gauge invariant information con- 
tained in the (untraced) Polyakov loop, which in lattice nota- 
tion is given by 



This is a temporal holonomy winding around the compact Eu- 
clidean time direction of extent N^. For gauge groups with 
nontrivial centre the traced Polyakov loop, 

L^=trF'Vx. (2) 

with the trace being taken in the defining representation^ 
serves as an order parameter for the deconfinement phase 
transition. The phase transition goes along with spontaneous 
breaking of the centre symmetry resulting from non-periodic 
gauge transformations under which 

L^^zL^, ze Z{Nc) . (3) 

The deconfined broken- symmetry phase at sufficiently large 
Wilson coupling /3 > /3c is characterised by (L) ^ (see llal 
for a review). 

Recently it has been found that gauge groups with trivial 
centre may also lead to a deconfinement transition depending 
on the size of the gauge group fT^flTlfTsll . 

At this point the choice of dynamical variables needs to 
be addressed. Under periodic gauge transformations g G 
SU {Nc) the holonomy transforms as 

'i^x^gx'i^xg^^ , (4) 

which leaves its eigenvalues and, as a consequence, its trace 
Q invariant. The eigenvalues are permuted arbitraril y by 
gauge transformations corresponding to Weyl reflections 11^ . 
This invariance is taken into account by constructing symmet- 
ric polynomials in the Nc eigenvalues. From unimodularity 



(1) 



We do not include a normalisation factor 1 /Nc for the ease of later nota- 
tion. 
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the product of the eigenvalues is 1 and there are only Nc — 1 
independent polynomials, for example the traces tr^?^^ for 
1 < n < Nc. These in turn are in one-to-one correspondence 
with the characters of the Nc — 1 fundamental representa- 
tions (see below). Hence, for SU {2) Lx (which is real) is 
sufficient while SU{3) requires Lx and L* , the latter being a 
linear combination of L^. and Ivf Only for Nc > 4 traces 
of higher powers of ^ are needed as independent dynamical 
variables 12011 . 

In two recent papers |El we have studied Polyakov 
loop models on the lattice for the simplest non-Abelian gauge 
group SU{2). The models have been derived using strong- 
coupling techniques at small Wilson coupling P and a newly 
developed inverse Monte Carlo (IMC) method which works 
for arbitrary values of (3. The latter method allows for a map- 
ping of Yang-Mills theory at a certain value of f3 to any appro- 
priately chosen Polyakov loop model of the form 

^efT= A/j(/5)(x/(^i^)xj(?/)+h.c.) , (5) 

{xy),IJ 

where the summation is over nearest neighbours and group 
representations /, J. The group character xi is the trace of 
the Polyakov loop in representation /, 

Xi{x) = xi['^x]=tri^a^ . (6) 

All characters xi polynomials in the characters corre- 
sponding to the fundamental representations. For SU{2) the 
simplest model is of Ising type. 

Si = X LxLy , (7) 

as first suggested in lEsIl (see also j^lz^]). 

The output of the IMC routines are the /^-dependent ef- 
fective couplings A/ J. Even without particular knowledge of 
these one may study the models ^ in their own right as statis- 
tical field theories. The Svetitsky-Yaffe conjecture may then 
be utilised to deduce information about the Yang-Mills phase 
transition. For SU (2) we have been able to show that the 
mean-field analysis of the effective models yields a surpris- 
ingly good agreement with the Monte Carlo analysis of both 
the model itself and the underlying Yang-Mills dynamics. It 
should be stressed at this point that the models (0 are not just 
simple scalar field theories with a linear target space as all 
characters ^ take values in a compact space parametrised by 
the Nc — 1 fundamental characters. Hence, in the lattice path 
integral each lattice site is endowed with the reduced Haar 
measure on the gauge group rather than a Lebesgue measure. 

This paper presents our first steps to generalise the results 
of (ill Q to the gauge group SU{?>). We label the char- 
acters by two integers p and q which count the numbers of 
fundamental and conjugate representations ('quarks' and 'an- 
tiquarks' in the 6'/7(3) flavour language) required to construct 
the representation (p, q). Equivalently, these integers charac- 
terise the horizontal extensions of 5'/7(3) Young tableaux (see 
Fig. [I}. Under the Z(3) centre transformations the characters 



P 



q ► 

FIG. 1: Character labels and 6't/(3) Young tableaux. 

transform according to the rule 

p *q p—q 

Xpq ~^ ^k^k Xpq = ^k Xpq i 

Zk=e^v{^k^ , fe = 0,l,2, (8) 

so that the most general centre- symmetric effective action 
with nearest-neighbour interaction may be written as 

^eff[x] = Y Kq.p'q' {xpq{^) Xp' q' {v) ^^'^^ 

{xy),pq,p'q' 
p-\-p' —q-\-q mods 

+ Y Apg,oo (xpg(a^) +h-c.) . (9) 

a3, pq 
p—q mod 3 

This coincides with the ansatz suggested by Dumitru et al. 
|[l4l] . The first sum in the effective action consists of hopping 
terms involving monomials of the form L^Ly or L^Ly^ (and 
h.c.) while the second sum is a 'potential' term containing 
only powers (and h.c.) localised at single sites. 

The remainder of the paper is organised as follows. In Sec- 
tioniniwe confirm the ansatz ^ by means of a strong coupling 
(small-/3) expansion for 5'eff[x]- For a restricted set of cou- 
plings (and hence representations) we investigate the resulting 
Z(3) models by minimising the classical action (Section Ullll 
followed by an improved mean field analysis in Section |IV| 
In agreement with the Svetitsky-Yaffe conjecture we find a 
first order phase transition from the symmetric to a ferromag- 
netic phase. Our improved mean field analysis already reveals 
an interesting phase structure with four different phases: two 
distinct ferromagnetic phases, one symmetric and one anti- 
ferromagnetic phase. Besides the first order transitions we 
detect second order transitions from the symmetric to an anti- 
ferromagnetic and to a ferromagnetic phase. The continuous 
transition from the symmetric to the anti-ferromagnetic phase 
is to be expected since the models reduce to the 3-state Potts 
model for Polyakov-loops having values in the centre of the 
gauge group. Section fVl contains the results of our extensive 
Monte Carlo simulations performed on a Linux cluster where 
we have implemented the powerful package jenLaTT. Simi- 
larly as for SU (2) the mean field and numerical results are in 
surprisingly good agreement. This is presumably due to the 
existence of (three) tricritical points. Depending on the or- 
der of the transition we localise the critical lines with either a 
Metropolis, a multicanonical or a modified cluster algorithm. 
In addition we have checked that the critical exponents u and 7 
at the second order transition to the anti-ferromagnetic phase 
agree with those of the 3-state Potts model in 3 dimensions. 
Finally, in Section|Vl|we wrap up with discussion and conclu- 
sions. 
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II. STRONG-COUPLING EXPANSION 

In this section we briefly recapitulate the strong-coupHng 
(smaU-/^) expansion for the SU (3) Wilson action at finite tem- 
perature . It is known that the leading order result 
(P^^ ) stems from ladder diagrams that wind around the tem- 
poral lattice extension and corresponds to an Ising type model 
analogous to 0. By going beyond the leading order we will 
encounter higher group representations/characters and hence 
have an independent confirmation of the ansatz ^ for the ef- 
fective action. 

Our starting point is the standard Wilson action, 



1 



Nc 



Re tr Ua 



(10) 



where the summation over plaquettes contains both temporal 
and spatial links. The effective action 5'eff[^] is introduced as 
usual by inserting an appropriate (group valued) delta function 
in the path integral, 



Z = / D/7e" 



(11) 



While it is not known how to perform the final integration an- 
alytically for the full Wilson action one can straightforwardly 
integrate order by order in p. Thus we expand the Boltzmann 
weight 



k 



(12) 



and integrate separately over temporal and spatial links, 
D/7 = TiUtTiUs. Adopting temporal gauge we set all tempo- 
ral links equal to unity apart from the links in the first timeslice 
which according to {I} may be identified with the Polyakov 
loop. Integrating out all spatial link variables we obtain the 
partition function 



Z = exp (^ogY^OkP^^ 



(13) 



where we have introduced the operators Ok = J DUsOk- The 
effective action may finally be written as 

5e,= -log(^0,^M =^^n/3^ (14) 
\ k In 

where the coefficients Sn are related to the operators via 
the linked-cluster theorem. In the remainder of this section 



we are going to determine the explicit form of the effective 
operators Sn . 

We first rewrite the Wilson-Boltzmann weight (IT2I1 as 



e '^^ = exp 



(15) 



and expand the single-plaquette contribution in terms of 
5(7(3) characters, 



(16) 



with / = (p, q) (see Fig.[l]). All P dependence now resides 
in the generalised Fourier coefficients aj which accordingly 
may be further expanded. 



(17) 



I,k 



An explicit computation of the shows that these vanish 
whenever the representation labels become sufficiently large, 
namely if I /| = p-\-q > k | 27]. This yields the important inter- 
mediate result that to any given order k in the strong-coupling 
expansion only a finite number of characters contributes, 

k \\I\=0 J 

The integrations over the spatial links are standard group in- 
tegrals which can be found in the texts 13 ES- The upshot is 
that only connected link arrangements ('polymers') wrapping 
around the temporal extent of the lattice yield nonvanishing 
contributions. The leading term is a ladder diagram consist- 
ing of Nr plaquettes each of which contributes a factor of /3 
implying a total contribution of 0{P^^). The associated op- 
erator is explicitly found to be 

On^ (X Xl0(?^a3)X0l(?^a3+^) +h.c. . (19) 

A typical operator of order P^^^ is given by 

E E^H/5)(x/(^^a3)xK?^a.+.)+h.C.) , (20) 

I,\I\=j x,i 

in terms of which the Wilson-Boltzmann weight JTSt becomes 



-Sw 



exp l^-lnX^On/?"^ 



(21) 



r=l.../c ai...ar 

o-iH \-ar-^k 



Expanding this to next- to-leading order (P'^^^ ) yields the sim- 
ple expression 

^eff = S^^S^^ S^S^ + 0{P^^^) . (22) 
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Note, however, that some care has to be taken in interpret- 
ing products such as S'^S^ which by (|2Q|) also contain discon- 
nected pieces. Upon expanding the logarithm by means of 
the linked-cluster theorem we are led to keep only connected 
contributions of the form 

S'S^^ E E(xH'Pa>)xK'Pt«+fc)+h.c.) 

I,\I\=i x,k 
J,\J\=j 

X (^XjiVx)x*jiVx+k)+h.c}j . (23) 
Making use of the character reduction formula 

Xi{x)xj{x) = ^ Cfj xk{x) , Cfj e R , (24) 

K 

products of characters at the same site may be reduced to sin- 
gle characters. As a consequence, the connected part of (l22l) 
takes the explicit form 

S., = Aio ^ (xi0(^a3)X0l(^a3+i) + h.C.) 

+ A20 E (x20('Pa=)X02('Pa=+i) + h.C.) 

X,i 

(25) 

+ A2l^(x20(^^a3)Xl0(^^a3+^) 

x,i 

+ Xw{'Vx)X2o{Vx+i)+h.c}j 
+ PiJ2xniVx)+0{/3'''^). 

X 

For what follows it is useful to introduce the short-hand nota- 
tion 

^eff = AioS'io+A2o5'2o+AiiS'ii+A2iS'2i+piVi+0(/3^^'^) , 

(26) 

with the obvious term-by-term identifications as compared 
to (l25t . Our conventions are such that all couplings in (l25t 
and (I2SI) are real functions of (3, the single leading one being 
Aio = 0{(3^^) (as noted already in |24]) while the sublead- 
ing ones are 0(/3^^^ ). It is straightforward to include higher- 
order terms the number of which increases rapidly. At order 
l^3N^ , for instance, there are already 1 1 terms so that we re- 
frain from going beyond next- to-leading order in (3. 

For later purposes it is useful to express the operators ap- 
pearing in (l26l in terms of the fundamental loops L and L*. 
Octet and sextet characters (xii and X20, respectively) are 
eliminated via the standard reduction identities 

3 (g) 3* = 8 © 1 and 3 3 = 6 © 3* , (27) 

which are equivalent to the character relations (recall that L = 

Xio) 

Xii = |Lp-l and X20 = i^'-i^*. (28) 



Making use of the latter the different terms in (l26l) become 

Sw = E i^^^l + h-c-) , (29) 

{xy) 

(30) 

5n = E {\L.?\Ly? - \L^? - \Ly? + l) , (31) 

S21 = E (LlLy - LlLy + LlL^ - L^L; + h.C.) , 

(xy) 

(32) 

V^ = ^{\L^\'-l) . (33) 

From these expressions it is obvious that each operator Spq is 
manifestly real. Under charge conjugation Xpg Xpq = Xqp^ 
hence Lx ^ L^, whereupon all terms in ^eff are charge conju- 
gation invariant as required (l^ . Note that the orders included 
correspond to terms that are quadratic, cubic and quartic in L 
and/or L*. Higher powers will arise upon taking into account 
higher representations. Thus, in this respect 6*2 1 is somewhat 
singled out being of only cubic order. This will become im- 
portant in a moment. 

III. QUALITATIVE CLASSICAL ANALYSIS 

It has already been pointed out by Svetitsky and Yaffe ^ 
that effective actions with Z(3) centre symmetry are closely 
related to the 3 -state Potts model which shows a first order 
phase transition from a symmetric to a ferromagnetic phase 
1301 131I [32I1 . To make the relation manifest we restrict the 
Polyakov loop to the centre elements Zk introduced in (|8|). 
Setting = we find the general formula 

Xpqi^k) = Zl'^^dpq , (34) 

where dpq denotes the dimension of the representation (p, g), 

= ^(p+l)(^ + l)(p + ^ + 2). (35) 

Applying this to the effective action in (I2SI) we find, up to an 
additive constant. 




kx e {0,1,2}, 



(36) 

with effective coupling 

A= 18(Aio + 4A2o+4A2i) . (37) 

The action is precisely the one of the 3-state Potts model 
mill. We thus expect that our effective Polyakov loop 
models will have a phase structure generalising the one of 
the 3-state Potts model. The latter is known to have a ferro- 
magnetic phase for large negative A and an anti-ferromagnetic 
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phase for large positive A |35|l. As a preparation for 
the discussion of later sections it is hence useful to obtain 
some qualitative understanding of the phase structure of the 
Polyakov loop models viewed as generalisations of the 3 -state 
Potts model. For the reason mentioned at the end of the previ- 
ous section we choose as a minimal generalisation the follow- 
ing effective action, 



10*->10 



21 



(38) 



which in terms of the fundamental loops may be written ex- 
plicitly as 



5e.= (Aio-2A2i)^ (L,L;+h.c.) 



+ A21 i^l^y + + h.c.) . (39) 

{xy) 

Note that there are also quadratic contributions stemming 
from 6*2 1. The action (l39l is manifestly Z(3) centre symmet- 
ric under Lx ZkLx- 

It is important to realise that (l39t differs from the standard 
lattice actions for scalar fields in several respects. First, the 
field Lx is dimensionless, being the trace of a unitary matrix. 
This allows for the presence of cubic hopping terms connect- 
ing neighbouring sites. Even more important is the fact that 
the target space of Lx is compact. Introducing the eigenvalues 
of ^ via 



(40) 



and writing L = Li + zL2we find for the real and imaginary 
part of L, 

Li = cos + cos 02 + cos(0i + ^2) , (41) 
L2 = sin + sin 02 - sin((/)i + ^2) • (42) 

The target space of L may then be sketched in the complex 
L-plane (see Fig.|2|). The boundary corresponds to the points 
with 01 =02, the singular 'corners' being given by the three 
centre elements ^ = Zkl. Let us try to get some first rough 
idea of the phase structure associated with the two-coupling 
model (l39l in the A10-A21 plane by looking at the extrema of 
the classical action. If we vary the couplings these will trace 
out a certain (possibly discontinuous) trajectory in the target 
space given by the triangle of Fig.|2| 

As we argued earlier, for centre- valued Polyakov loops the 
effective action (l39l) reduces to the action of the Potts model 
(l36l with coupling A = 18(Aio + 4A21). Thus we expect 
a ferromagnetic phase (F) for large negative Aio + 4A21 and 
an anti-ferromagnetic phase (AF) for Aio + 4A21 large and 
positive. In a region around the origin in the coupling plane 
entropy dominates energy and we cannot expect to actually 
obtain the correct phase- structure in this region by purely clas- 
sical reasoning based on minimising the energy. Qualitatively 
we expect a symmetric phase in a neighbourhood of the ori- 
gin. This is represented schematically in Fig.|3lby the central 
rectangle. 



Im 




FIG. 2: Target space of the Polyakov loop L in the complex L-plane. 
The corners represent the three centre elements. The intermediate 
points (denoted anti-centre elements) will also become relevant for 
the discussion of the phase structure. 



In order to study the ordered phases (in particular AF) 
we divide the lattice in two sub-lattices (denoted 'even' and 
'odd') where the Polyakov loop takes values L^ and Lq, re- 
spectively. Two nearest neighbours belong to different sub- 
lattices. The absolute minima of the classical action 



5'eff(Le,Lo) oc (Aio - 2X2i){L^ 

+ A2l(L^Le 



L:+h.c.) 

^ LILo + h.c.) 



(43) 



will then be located at certain values Lq and Zq of the 
Polyakov loop which are identified with its 'expectation val- 
ues'. We trust this reasoning as long as we are sufficiently far 
from the origin of the coupling plane i.e. from the disordered, 
entropy-dominated phase. 

Any ferromagnetic ordering will be characterised by a min- 
imum with Lq = Lo = L while in an anti-ferromagnetic 
phase Lq ^ L^. Quite interestingly we find two distinct fer- 
romagnetic phases, one for which the Polyakov loop is near a 
centre element or L in the vicinity of Sz/c and a different fer- 
romagnetic phase with L taking values near the intermediate 
points marked by triangles in Fig. |2l We call this an anti- 
centre phase (AC). We expect a phase transition line separat- 
ing the ferromagnetic and anti-ferromagnetic phases at van- 
ishing Potts-coupling Aio + 4A2i. The resulting qualitative 
phase diagram is depicted in Fig.[3l 

To discuss the ferromagnetic phases it suffices to minimize 
the action (l43b with L^ = Lq = L,m which case 

S,,{L) (X (Aio - 2A2i)|Lp + 2A2i(L^ + L*^) . (44) 

This can be done analytically. To localise the anti- 
ferromagnetic phase we have calculated the absolute minima 
of (l43t on the target space depicted in Fig.|2|numerically. The 
combined analytical and numerical results are summarised as 
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-0.4 




FIG. 3: Qualitative prediction of the phase diagram in the couphng 
constant plane for the effective action (l39l . The ferromagnetic (F), 
anti-centre (AC) and anti-ferromagnetic (AF) phases are obtained by 
looking for classical minima. The symmetric, disordered phase (S) 
is located where entropy is expected to dominate over energy. 



follows. For negative A21 we have a transition 
while for positive A21 there is a richer structure, 



(45) 



p Aio^-3.1962A..^ Aio— ttA.i^ ^p ^^^^ 

The behaviour of a suitably projected order parameter £r (the 




FIG. 4: Behaviour of the order parameter ir defined in (l97t as a 
function of Aio for fixed A21 = 1. For comparison we have added 
the result from the mean field analysis to be developed in the next 
section. 

precise definition of which will only be needed later on) for 



positive A21 = 1 is shown in Fig. |4l Upon inspection one 
notes that for Aio sufficiently negative the system starts out 
with the Polyakov loop at a centre element. Increasing Aio 
beyond — 7A21 the order parameter drops monotonically until, 
at a critical coupling Aio ^ — 3.1962A21, there is a jump to 
the AC phase with L near a anti-centre element. The jump of 
ir is due to a centre-transformation and does not imply that 
the Polyakov loop itself jumps. Indeed, L changes smoothly 
and arrives at an anti-center element for Aio = — 3A21. The 
system stays there until Aio = — 28A21/II, where it jumps 
again, this time to the AF phase. As expected, we see no 
symmetric phase in a purely classical analysis. Actually, for 
A21 = 1 there is just no symmetric phase. 



IV. MEAN FIELD APPROXIMATION 

The next step of refinement to be presented in this section 
is a mean field (MF) analysis of the effective action (l39t . This 
will serve as a basis for a comparison with results from di- 
rect Monte Carlo simulations to be discussed later on. Due to 
the peculiarities of the model as compared to standard scalar 
field theories the application of the MF approximation is not 
entirely straightforward. For the benefit of the reader we will 
set the stage by giving a brief outline of the necessary modi- 
fications. For further details the reader is referred to our ear- 
lier paper 1I22I1 . To keep the discussion sufficiently general 
we will first treat the effective action (l26l) with five couplings 
focussing on simpler examples later on. 

We are interested in expectation values which are computed 
by evaluating integrals of the form 



(A) 



ml 



Am 



(47) 



which apparently extend over the whole group manifold em- 
ploying the Haar measure dji^^x)- However, due to the gauge 
invariance of both action and measure the integrals can be 
reduced to the coset space of conjugacy classes which we 
(somewhat symbolically) denote by Vx • Hence we integrate 
with the reduced Haar measure by replacing 



Thus, (l47t is equivalent to 



(A) 



VP = l[d^red{-P^) ■ 



(48) 



(49) 



The probability measure W exp(— 5'eff)/Z^[0] is character- 
ized as the unique solution to the variational problem 



inf (S'eff + logp)p 
p 



logZ[0] = -T^[0] . (50) 



The expectation value on the left-hand side has to be taken 
with respect to the integration measure JJ^ dfiy^edCPx) P^P] 
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with p[P] denoting the probabiHty density of V. From this 
point of view the MF approximation is nothing else but the 
restriction of the permissible densities to product form, 



(51) 



Expectation values can now be simply computed site by site 
via factorisation, 

{Xl{r.)Xj{Vy)) ^ {Xl{ra:))Mj{'Py))y ■ (52) 

Our goal is to compute the effective potential as function of 
the mean characters Uj^f = ^mf(x/)- For that purpose we 
solve the variational problem (l5Ql on the space of product 
measures with fixed expectation values of the characters. This 
is done by introducing appropriate Lagrange multipliers jj. 
For ferromagnetic systems one may assume that the weight 
functions Px at each site are identical, Px = P- This assump- 
tion corresponds to a translationally invariant ground state. 
According to the discussion of the previous section we ex- 
pect anti-ferromagnetic phases and hence we must refine our 
choice for px . We therefore introduce different weight func- 
tions on the even and odd sub-lattices, respectively, 

IPo(P.) :sgn(x)=-l, ^^^^ 



defining the sign of a lattice point as 

sgn(a^) = (-1) 



(54) 



As a consequence, expectation values of characters will sub- 
sequently have two values depending on the sub-lattice where 
they are evaluated, 

[X/,o : sgn(a?) = -1 . 

The sources are taken constant as well when restricted to the 
even and odd sub-lattices, jj^x = ji,e or jj^o, respectively. 
Like the characters the sources are complex. 

The action (l2Tl couples only nearest-neighbour sites so that 
its expectation value entering (l50ll may be written as 

(5'eff) = Vd |Aio (xio,oXoi,e + h.c.) H | 

V _ 

+ -^Pi (Xll,o + Xll,e) , (56) 

with V = denoting the lattice volume in d spatial dimen- 
sions. The logarithm in (BOll decomposes as 

V 

(logp) = — ((logpo)o + (logPe)e) • (57) 

It is convenient to drop the common volume factor V and con- 
sider densities instead. The variation of (ISOb finally yields the 
weight function for the even sub-lattice. 



Peir) 



1 



exp{-piVi{r) 

-j,-x{V)+jl-X*{V)} . (58) 



and a completely analogous expression for the odd sub-lattice. 
Here we have introduced j • % as a short-hand for ji xi 
and the single- site partition function 

f) = J d^rediV) exp { - piVi{V) 

+ . (59) 

The sources ^ are eliminated by inverting the relations 

d * _ * 9 

(60) 

to be evaluated separately on both sublattices. The Sch winger 
function w{j^j*) is defined as usual, 

= iogz(i,i*) . (61) 

Introducing the Legendre transform of (l6Tl according to 

7o(%, %*) = mf {j • X + j* • %* - (62) 

the solution of dSOl is finally obtained as the MF potential 
(density) as a function of even and odd mean fields, 

Um^{Xe^ Xt^ Xo^ %o) = d |Aio(xio,oXoi,e + h.c.) + . . . I 

+ ^Pl(Xll,o + Xll,e) + ^7o(%o. %o) + ^7o(%e. %e) • 

(63) 

From this expression one can easily derive relations between 
the sources j and mean characters %. For instance, by varying 
^mf with respect to x 10,0 we obtain 

= d (AioX01,e + A2iX20,e) + ^jio,o , (64) 

where we have used that the current is given as 

d 



JlO,o 



SXlO,o 



7o(Xo.%o) • 



(65) 



The first term in (l63t derives directly from the effective action 
(I2SI and hence contains four couplings to which the potential 
coupling pi has to be added. A complete MF analysis of this 
system becomes very awkward. In what follows, we therefore 
specialise to effective actions with only one or two couplings. 



A. One Coupling 

The action (l29b defines what we call the 'minimal model' 
with one coupling only, 

= XSio = A ^ {Lo^l; + h.c.) . (66) 



8 



Identifying (xio)e,o = ^e,o and inserting (l66l) the MF poten- 
tial (I53l simplifies to 

ix^f (Le, L:, Lo, = d A (LeL: + LoLl) 

+ i70(i^e,i^e*) + ^70(i^o,i^:). (67) 

It is useful to define the generic MF potential 

v^,{L, L'') = d\ |Lp + 7o(L, L*) , (68) 
so that (l67l can be rewritten as 

'^mf(^e, Ll.Lo, L*) = - — |Lo - Le|^ 

+ \^,{L,, LI) + i^^f (Lo, LI) . (69) 

This expression clearly shows that for negative A configu- 
rations with Lo = Lq are favoured making the distinction 
between even and odd sub-lattices obsolete. The remaining 
unique expectation value, 



For arbitrary order parameters we have solved the gap equa- 
tions (l73l numerically. The result is depicted in the follow- 
ing figure. The transition S-AF is second order, the one to 
the phase F first order. Thus, the anti-ferromagnetic transition 
must be at A+ = l/2d. The first order ferromagnetic transi- 
tion is slightly above — l/2d. Our MF analysis thus confirms 



3.0 



2.5 



2.0 



(70) 
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thus serves as an order parameter for the ferromagnetic phase 
transition. On the other hand, if A > Lg and Lq will cease to 
be positively correlated and we expect an anti-ferromagnetic 
phase transition. In this case, a reasonable order parameter is 
given by (ISIll 



Let us also introduce the absolute values of the order parame- 
ters, henceforth denoted by 



I = \L\ and m = \M\ . 



(72) 



If we assume for the moment that the occurring phase tran- 
sitions are second order the corresponding MF critical cou- 
plings can be computed analytically as we are going to 
demonstrate next. The consistency condition (IS4l reduces for 
the case at hand to 



= dXLo + ^ je and = dXL^ + ^ jo 



(73) 



As both Lq and Lq may be assumed small near the criti- 
cal coupling Ac the corresponding sources will also be small. 
Thus, we may expand 

^^ /^M.e.exp{y,o+hx.}xio ^^.^ . (74) 
J dfired explj xio + h.c.j 

which again holds separately on each sublattice, Lq and 
I/O ^ jo- Plugging this into (l73t yields {2dX)'^ = 1 and hence 
the critical couplings 



A+ = ±. 



'2d 



(75) 



FIG. 5: Mean field results for the minimal model with one coupling. 
The first order ferromagnetic transition S-F is at A- = —0.13433 
and the second order transition S-AF at A+ = 0.166667. 

the qualitative results from the preceding section that already 
in the simplest model there is both a ferromagnetic and an 
anti-ferromagnetic transition. This is qualitatively consistent 
with the phase diagram Fig.jSjrestricted to the horizontal axis. 



B. Two Couplings 

Even for the simple model of the previous subsection no ex- 
plicit expression is known for zq {j , i* ) as the SU (3) group in- 
tegrals cannot be evaluated in closed form (a fact well known 
from strong-coupling expansions, see e.g. 12^ ). Things natu- 
rally become worse if additional couplings are turned on. Al- 
ready for two couplings i.e. the action (l39t , the only way to 
proceed is by means of numerical methods. 

In order to obtain the MF version of the phase diagram 
Fig. [3] we have employed the following algorithm: 

1. At the extremal points of (IS3t all sources jo,e occurring 
in (l58b can be eliminated in favour of the expectation 
values Xo,e as in (IS4l) . Since the character target space 
is compact it can be easily discretised defining measures 
Pe and po at each point. Using these measures expecta- 
tion values (x)o,e can now be computed which in gen- 
eral will differ from Xo e- Hence, we first look for local 
minima of 

'^(Xo,Xe) = ll(x)o-Xoll + ll(x)e-Xell, (76) 

with norm \\v\\ = \vi\. These minima, however, do 
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not correspond to exact solutions yet but rather serve as 
the starting points of a recursion. 

2. We now solve the equation cr{Xo^Xe) = for and 
Xe by Newton iteration using the local minima of the 
previous step as initial input. In this way we end up with 
multiple solution vectors Xo e ^^^^ extremising (l63l . 

3. The solution vector with minimal Umf contains the de- 
sired ground state expectation values. 

With the help of this algorithm we are able to compute the 
expectation values of L (TTQI and M (1771 . 

We conclude this section with some remarks concerning 
our choice of the sources j. Allowing them to be complex 
yields an 8-dimensional parameter space for the observables 
of model (l39l with two characters, xio and X2i- For this 
large parameter space it would actually be simpler to perform 
a high-precision Monte Carlo simulation than to find a good 
approximation for the global minima of (l63l . For this reason, 
we have chosen real sources as an input to our MF approxi- 
mation. We expect this to be a very good approximation as 
long as the peak of the probability distribution for xio is con- 
centrated near the real axis in Fig.|2| 

Comparing with the classical analysis we note that our real- 
source assumption is justified for all couplings which are lo- 
cated away from the boundary between the phases F and AC. 
The F-AC transition should be second order according to the 
analysis of Section Hill while the MF approximation with real 
sources predicts a first order transition (see Fig.|4|l. The con- 
tradiction will be finally resolved by Monte Carlo simulations 
showing that near the F-AC transition the peaks of the prob- 
ability distribution for xio are off the real axis. Apart from 
this fairly small region in parameter space we find a remark- 
able agreement between the MF results and the Monte Carlo 
results of the next section. This provides the ultimate justi- 
fication for our simplifying choice of real sources in the MF 
approximation. 



represent the conjugacy class according to 

m ^ V.i^) = diag(e^^Se^^%e-^(^^+^^)) , ^^^^ 

Here, the following restrictions should be imposed such that 
the angular coordinates cover each class only once, 

0< 01,02, 01 <02, 02 < (-0i-02)mod27r. (78) 

These restrictions are somewhat awkward to implement in a 
simulation code. It is much more convenient to let $ take 
values in the full square [0, 27r) x [0, 27r) which covers the 
fundamental domain given by (l78l six times. Due to the resid- 
ual gauge symmetry of the system it is clear that expectation 
values will be unaffected by this over-counting. 

In the coordinates (1771 the reduced Haar measure becomes 

^ fif.\ ^ . 2 01 - 02 A 

X sin^ (^01 + sin^ (^02 + d0i#2 , (79) 

where the normalization is such that the measure integrates 
to unity over the square [0, 27r) x [0, 27r). All characters can 
be expressed in terms of 0i and 02 . Using (I79l the full mea- 
sure (l49l may be straightforwardly expressed in terms of the 
angular coordinates as 

VV e-^^^^[^] = II d/ired{^x)e-^^''^^^ . (80) 

X 

The focus of our numerical studies has been the phase dia- 
gram in the Aio — A21 plane. There, we have scanned through 
the region [-0.25, 0.33] x [-0.22 . . .0.16] with a resolution 
of 71 X 46 points, which were in total 3266 different Monte 
Carlo simulations. Before we present our results a few words 
on our numerical techniques are in order. 



V. MONTE CARLO SIMULATION 

The previous two sections have provided us with a good 
deal of information on the phase structure of the effective 
model in the most interesting regions of parameter space. 
Based on this we have performed a large number of Monte 
Carlo simulations to quantitatively check the MF predictions 
and obtain a precise picture of the critical behaviour. As be- 
fore, to avoid excessive complexity, we have concentrated on 
the action (l39l with two couplings Aio and A21. 

At the beginning of Section |IV| we have already noted that 
both the action and the reduced Haar measure depend only on 
the conjugacy class of the (untraced) Polyakov loop. Hence, 
one can choose the basic field variables either as the traces in 
the fundamental representations (L, L*) and powers thereof 
or as a suitable parametrisation in terms of the eigenvalues as 
introduced in (l4Ql . It turns out that for a numerical treatment 
the latter proves to be more appropriate and so we use (l4Ql to 



A. Algorithms 

For the investigation of phase transitions, in particular their 
order, histogram methods are widely used and accepted. This 
approach, however, requires large statistics and thus tends to 
consume a lot of computer time. In addition, we are inter- 
ested in a fairly large range of coupling constants. For these 
reasons, the updating algorithm for the Polyakov loop mod- 
els has to be fast and versatile. It turns out that the standard 
Metropolis algorithm favourably matches both requirements 
if we aim at an accuracy of about 5% — 10%. On the other 
hand, because of the highly nontrivial probability measure in- 
volved, a heat bath algorithm does not seem applicable or, in 
any case, would be too time consuming. In addition, its local 
nature should not yield any enhancement of statistics near a 
first order phase transition. We have thus refrained from im- 
plementing a heat bath update scheme but rather decided to 
optimise the Metropolis algorithm as described in the follow- 
ing two paragraphs. 
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FIG. 6: Histogram of i with Aio = -0.13721 on a 10^ -lattice sam- 
pled against the histogram for = 9. 



1. Multicanonical algorithm 

When a system undergoes a first order phase transition, the 
histogram associated with the order parameter will generically 
display a multi-peak structure. Depending on the total volume 
of the system the peaks can be very pronounced. In other 
words, the configuration space is decomposed into distinct 
sectors between which local algorithms can hardly mediate. 
One way to overcome the resultant failure in sampling the to- 
tal configuration space is to make use of the multicanonical 
algorithm, see e.g. f^]. 

The crucial improvement step consists in the replacement 
of the measure used in (ISOl) according to 



d fired exp(-S') djired exp(-S') r]{l) 



(81) 



The new improved measure on the right-hand side has a 
weight T] = r]{i) which depends on the (modulus of the) order 
parameter. One chooses the particular form 



(82) 



where p{t) denotes the probability density of the order pa- 
rameter. This choice leads to an enhancement of configura- 
tions that would otherwise be suppressed and thus allows for 
a much improved ergodic behaviour of the algorithm. 

The effect due to the altered measure is illustrated in Fig.|6| 
where two typical distributions are plotted. In the original 
distribution one clearly recognises two well- separated peaks. 
Any local algorithm will fail to sample such a distribution 
properly once it is trapped at one of its peaks. The distri- 
bution actually used 'closes the gap' enabling transitions be- 
tween different peak regions during the simulation. In the end, 
of course, one has to correct for the change in the measure by 
reweighting with 7/~^, 



{Qri-\i)) 



mult. 



mult. 



where we have denoted expectation values taken with respect 
to the modified measure with a subscript 'mult' . 

A slight problem with this approach, however, still has to 
be overcome. One actually needs right at the beginning what 
one set out to compute originally, namely the distribution p. 
A practicable strategy is e.g. the following. From a small lat- 
tice volume, say Vq, where peaks are usually less pronounced, 
one obtains an approximate distribution function po (^) which 
is then used on a slightly bigger lattice, say of volume Vi. 
This simple trick can be further refined if on the larger lattice 
one first computes pi {£) using po and subsequently repeats 
all measurements employing pi . In practice, this procedure is 
iterated several times to make larger and larger lattices avail- 
able. Going beyond a volume of F = 12^ requires additional 
knowledge of the scaling behaviour of p{£). Fig.Qshows that, 
to a very good approximation, the scaling depends linearly on 
the volume, V = N^, 

log p{i, N) ^ A{i) + C{i) . (84) 

In summary, the multicanonical algorithm yields substantial 




(pie)) 



mult. 



(83) 
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FIG. 7: Plot of log p{i, N) against volume V = N for the minimal 
model with Aio = —0.13721 confirming the scaling behaviour i84t . 
The coupling is very close to critical. 

improvements compared to the standard Metropolis algorithm 
and allows for very accurate simulations. Lattices with vol- 
umes up to F = 20^ could thus be studied near the first order 
phase transition. However, as the implementation is very in- 
volved we applied it only to the minimal model (A21 = 0). In 
principle, the generalisation to incorporating further couplings 
is straightforward, being merely a matter of having sufficient 
computing time available. 



2, Cluster algorithm 

The multicanonical algorithm is best suited for studying 
first order transitions. This is no longer true for second or- 
der transitions where it is outperformed by cluster algorithms 
- at least if there is one available. For our purposes, none 
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of the algorithms on the market could be immediately put to 
use. We therefore decided to modify the well-known Wolff 
algorithm |3^, an extension of the Swendsen-Wang Is^ al- 
gorithm originally proposed for discrete spin systems. 

Before actually describing our modifications let us briefly 
recall the main idea of Issll . Suppose the action (and the mea- 
sure) are invariant under a certain global symmetry which acts 
on the fields according to 



(85) 



Typically, the symmetry operator R will depend on some pa- 
rameters which we collectively denote by co (continuous or 
discrete) so that R = R{uj). It is important to note that the 
operator R has to be idempotent, R'^ = id, in order to ensure 
detailed balance. The algorithm then works according to the 
following list. 

1 . Fix some parameter ujq and hence some transformation 
Rq = R{coo). Randomly choose a lattice point x and 
apply the symmetry transformation 



= Ro 



(86) 



which may be viewed as flipping the field variable 
The point x is checked and added to the cluster. 

2. Repeat the following for all unchecked neighbours y of 
x: 

Let S = S{^x^^y) denote the contribution to the ac- 
tion from the link {xy) and compute 

= S{^'^, Ro ^y) - Si^'x. ^y) , (87) 
which may be rewritten as 

AS = S{^a.,^y) - S{^'^,^y) , (88) 

since S is already invariant under R. The decision to 
add y to the cluster is subject to an accept/reject step so 
that the probability to flip becomes 



AS>0, 
A5'<0. 



(89) 



If is flipped, check the point y. 

3. Go back to Step Q for all sites added to the cluster in 
the previous step. 

Note that ^ measures whether it is advantageous to flip 
once ^x has been flipped. From the very construction of the 
algorithm it should be already clear that the clusters will in- 
crease with the correlation length of the system. In this way 
one suppresses the phenomenon of critical slowing down ob- 
served with local algorithms. This makes the cluster algorithm 
particularly suited for the study of second order phase transi- 
tions. 

The first step in adapting the cluster algorithm to our needs 
is to find suitable symmetries of the action Sill . Being a 



Polyakov loop model the symmetry in question is the discrete 
Z(3) symmetry, L ZkL. In addition, the action has to be 
real which implies constraints on the way complex conjuga- 
tion acts. From these symmetries one can construct three op- 
erators i = 0, 1, 2, acting on the Polyakov loop according 
to 



RiL= {ziLf , Zi e Z(3) 



(90) 



As required the Ri square to unity. Furthermore, it is easy 
to see that both the operators Spq appearing in (l26l and the 
domain of the Polyakov loop L are left invariant by the action 
of Ri. The latter is illustrated in Fig. [Sl for a particular value 
of L. For the actual algorithm the transformations (l9Ql are not 
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FIG. 8: Illustration of the Z(3) reflections Ri used in our modified 
cluster algorithm. 

immediately applicable since the simulation is based on the 
angular variables $ introduced in (TTTl . However, it is just a 
matter of a little algebra to show that the Ri act on the angles 
^ via 



2) 



(27r — ^1, 27r — ^2) under Rq , 

- 01, ^ - ^2) mod 27r under Ri , 

- 01 , ^ - ^2) mod 27r under R2 . 

(91) 

Whereas in the original cluster algorithm flipping along ran- 
domly chosen lattice sites is sufficient to guarantee ergodicity 
this is no longer true in the present case. Hence, we have to 
augment the update scheme by standard Metropolis sweeps to 
make the algorithm ergodic. For one Monte Carlo step our 
cluster algorithm finally can be summarised as follows: 

1 . Choose a random number Nm between and V = . 

2. Do Nm standard Metropolis sweeps at randomly drawn 
lattice points. 

3. For a suitable fixed number A^ci repeat the steps for 
building a cluster as described above. 
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4. Do V — Nm additional Metropolis sweeps, again at ran- 
domly chosen lattice sites. 

From several test runs we have found that the number N^i 
should be chosen such that total number of flipped sites after 
performing step (3) is approximately half the number of all 
lattice sites. Thus, if |C| denotes the typical size of a cluster, 
the following equation should hold (at least approximately). 



V 

2\C\' 



(92) 



One of the most interesting questions of course is the gain in 
performance compared to e.g. the standard Metropolis algo- 
rithm. We have found that the autocorrelation time for the 
order parameter r£ is independent of both the lattice extent 
(at least in the range N = 8 ... 28) and the coupling con- 
stant (if close to the second order phase transition) implying 
a dynamical critical exponent of z = 0. Moreover, for an op- 
timal choice of Nc\ the autocorrelation time is of order unity. 
As a result, our cluster algorithm outperforms the Metropolis 
algorithm even for small lattices. On the largest lattices we 
have considered the cluster algorithm reduces autocorrelation 
times by two orders of magnitude as compared to Metropolis 
updating. This improvement comes at the cost of a slightly in- 
creased complexity, specifically a factor of 1.5 in computing 
time which clearly is negligible. On the other hand, in line 
with our expectations, no significant improvement has been 
found near the first order transition. 



B. Results 

We present the results of our Monte Carlo simulations in 
the same two steps as for the MF approximation. Hence, we 
first report on the minimal model and switch on A21 later 
on. With only A 10 different from zero it is reasonably cheap to 
perform highly accurate measurements so that a precise quan- 
titative comparison with the 3-state Potts model is possible. 

For the model (l39t with two couplings we determine the 
phase diagram in the A10-A21 plane and compare with our 
expectations as laid out in the previous two sections. We con- 
clude with a careful study of the nature of the phase transi- 
tions, in particular their continuity properties. 
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FIG. 9: Phase transition of Fig. fTol when calculated within MF ap- 
proximation. 
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FIG. 10: Ferromagnetic phase transition computed with the multi- 
canonical algorithm on a 16^ -lattice. The expectation value of i is 
plotted against its probability distribution given by the area shaded 
in grey. The latter clearly shows the correct discontinuous behaviour. 
Note that, since we measure the modulus, statistical fluctuations 
manifest themselves in a (small) positive value of the order parameter 
even in the symmetric phase. 



1. Minimal model ( one coupling ) 

To determine the ferromagnetic phase transition for Aio < 
we have used standard techniques which need no further ex- 
planation. It suffices to note that in order to study the first or- 
der (S-F) transition we performed 10^ sweeps on a 16^-lattice 
in the multicanonical ensemble for each value of the coupling 
constant. This led to highly accurate statistics until we ap- 
proached the close vicinity of the critical coupling itself. In 
this regime our statistics is restricted to 10^ independent sam- 
ples due to large correlation times of the order of 10^ sweeps. 
As Fig. [10| shows this is sufficient to demonstrate that the tran- 
sition is first order. 



In Fig.|9lwe compare our Monte Carlo result for ^(Aio) with 
its MF approximation. The figure basically zooms into that 
part of Fig. [S] where the first order ferromagnetic transition 
is located. Again, the first order nature of the transition is 
corroborated. 

In addition to the coupling dependence of the order param- 
eter we have determined both the critical coupling and the 
discontinuity A£ at Aio,crit- The results are given in Table |I| 
together with the MF prediction. Again we find a surpris- 
ingly good quantitative agreement between simulations and 
MF approximation. Moreover, if we consider the ratio of the 
two critical couplings, say Aio,crit;F/Aio,crit;AF, we are able to 
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compare this with results of the 3 -state Potts model. The ac- 
tual figures turn out to be fairly close, namely —0.6904 for 
the minimally coupled Polyakov model and —0.6750 for the 
3-state Potts model HEM- 



TABLE I: Critical couplings for the S-F and S-AF phase transitions 
and jump A^. For the first order transition (S-F) there is excellent 
agreement between MF and Monte Carlo data. Even the values for 
A£ agree within 10%. For the second order transition (S-AF) the 
critical couplings agree within approximately 20%. 



method 


AlO,crit (S-F) 


Ai 


AlO,crit (S-AF) 


Monte Carlo 


-0.13721(5) 


1.33(2) 


0.19875(5) 


Mean Field 


-0.13433(1) 


1.46(1) 


0.16667(1) 
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FIG. 1 1 : Expectation value of m near the anti-ferromagnetic phase 
transition computed via MF approximation. 

Table m also displays the (positive) critical coupling for the 
second order AF transition. This has been analysed with our 
modified cluster algorithm by performing a total number of 
2 X 10^ sweeps. With these large statistics at hand it is also 
possible to determine some of the critical exponents thus prob- 
ing the universality properties of the model. To do so we have 
employed standard renormalisation group techniques follow- 
ing 1401] . In particular, we consider the Binder cumulant U 
and susceptibility x given by 



U = l 



(93) 
(94) 



with m as defined in (ItTI and N denoting the spatial extent of 
the lattice as before. 

The Binder cumulant U = U {N^ Aio) is constructed such 
that it becomes independent of TV close to the critical point. 
Hence, the latter is rather precisely determined as the point 
where the graphs of U (plotted for different N) intersect. This 
behaviour is nicely exhibited in Fig|T3l 
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FIG. 12: Expectation value and probability distribution of m near 
the anti-ferromagnetic phase transition obtained from Monte Carlo 
simulations. To identify clear signals we have chosen a large lat- 
tice with V = 28^ and evaluated 5 x 10^ sweeps. In contrast to 
Fig-Golno discontinuity is observed. Again the expectation value of 
the symmetric phase is biased since we measure the modulus of the 
order parameter. 
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FIG. 13: Binder cumulants U for different lattice sizes as a function 
of the coupling Aio- The critical coupling Aio,crit is determined via 
the intersection point. 



From the standard relations at criticality, 

x(Aio,crit) ^ N^^^ , 



dU{N,X,o) 



dX 



10 



(95) 
(96) 



we have finally computed the critical exponents 7 and u which 
are listed in Table ITTI 

The uncertainty in Aio,crit quoted in Table |I| is mainly due 
to the fact that the different cumulants do not precisely meet 
in a single intersection point (cf. again Fig. [T3t . The error 
in the critical exponents is estimated from a least square fit 
to the logarithm of (l95t and (l96l) . To cross-check our results 
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TABLE II: Critical exponents for the second order AF transition of 
the minimal model. 



exponent 


3-state Potts [351 


minimal Polyakov 




0.664(4) 


0.68(2) 


7/zy 


1.973(9) 


1.96(2) 



we have measured the expectation value of m and its prob- 
ability distribution p{m) in analogy with the ferromagnetic 
transition already discussed. As in the former case, the ex- 
pectation value of the order parameter alone does not suffice 
to decide on the order of the transition. However, the proba- 
bility distribution shown in Fig.[T2|is quite different from the 
one in Fig.[lO| No discontinuous behaviour is observed now 
which provides further (numerical) evidence for a second or- 
der phase transition. Equally important, the critical coupling 
obtained is compatible with the results presented in Table |l| 
For the sake of completeness Fig.[n]shows the MF prediction 
for the transition S-AF. 

Comparing with earlier results on the 3-state Potts model 
we draw the important conclusion that the minimal 
Polyakov loop model is in the same universality class. 



2. The phase diagram for two couplings 
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FIG. 14: Fundamental domain T of the order parameter L obtained 
by identifying Z(3) copies according to the depicted arrows. 

Having discussed the minimal model at length let us con- 
tinue by switching on the second coupling A21 in order to anal- 
yse the phase diagram in the coupling constant plane. This re- 
quires a suitably chosen 'indicator' to distinguish the (at least) 
four phases (S, F, AC and AF) we expect in accordance with 
our MF analysis of Section |IV| While I and m clearly are 
order parameters for the minimal model they are numerically 
less suited for the model with two couplings. It turns out ad- 
vantageous to construct a new observable denoted which 



may be obtained from I by the following procedure. We first 
divide the domain of L into six distinct parts as shown in 
Fig -El The light- shaded region represents the preferred locus 
of the Polyakov loop in the ferromagnetic phase F, whereas the 
dark- shaded region corresponds to the anti-centre ferromag- 
netic phase AC. To eliminate the (numerically superfluous) 
Z(3) symmetry the first step in our projection is to identify 
the regions as indicated by the arrows in Fig.[T4| In this way 
we end up with a fundamental domain for the Z(3) sym- 
metry centred along the real axis. Every L is mapped into 
by a centre transformation. To finally obtain a real observ- 
able we project the transformed L onto the real axis. This 
projection results in the variable Ir the sign of which clearly 
distinguishes between the two ferromagnetic phases, < 
indicates the AC phase, > the ferromagnetic phase F, 
while £ = in the symmetric phase S. Mathematically, the 
projection of L to is given by 



Ir = < 



ReL 
-iReL 
iReL 



^ImL 



^ImL 



(97) 



To detect the AF phase we simultaneously measure m so that 
we can finally discriminate between all possible phases. 

Our main results for the phase diagram are shown in Fig.fTSl 
obtained by the modified MF approximation explained earlier 
and Fig. [T6l which displays the Monte Carlo data and hence 
constitutes the most faithful representation of the phase struc- 
ture. 

The Monte Carlo simulations were carried out on a 8^- 
lattice. The number of sweeps was chosen such as to reduce 
the jackknife error in the estimate of £r below 0.1. The in- 
dependence of our samples was ensured by demanding that 
the autocorrelation time r^^ associated with the observable £r 
was less then one percent of the total number of sweeps. As 
a result our simulations included at least 4 x 10^ sweeps far 
away from the critical regions and more than 10^ sweeps in 
their vicinity. 

It is reassuring to note that the qualitative phase diagram 
Fig. predicted from energy-entropy arguments is quantita- 
tively confirmed by both Figs. [21 and [T6| Comparing the 
latter in some detail it is once more remarkable how good the 
MF approximation works. Within large regions of parameter 
space it agrees with the 'real' data within 10% or less. Inter- 
estingly, the observable £r also seems to be sensitive to the AF 
phase (see lower right part of Figs. (TSlandfT^. However, any 
further discrimination between phases F and AF by means of 
£r is impossible. To lift this degeneracy one clearly needs the 
AF order parameter m as an additional input. 

It remains to be discussed why £r is sensitive at all to AF or- 
dering. In what follows we will provide a heuristic answer in 
the context of the 3-state Potts model. By definition, all spins 
are (anti)aligned in the (anti)ferromagnetic ground state. In 
the Z(2) symmetric Ising model with only two spin states each 
possible ground state is two-fold degenerate, independent of 
the particular ordering (F or AF). The counting of degenera- 
cies, however, is totally different for systems with more spin 
states (hence higher symmetry). In the AF phase the addi- 
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FIG. 15: Phase diagram of the model with two couphngs as obtained via MF approximation. The figure shows a contour plot of the ground 
state expectation value of ^r in the A10-A21 plane. The phase transition between the two ferromagnetic phases (F-AC) is clearly visible in the 
lower left part. Note that even discriminates between symmetric and anti-ferromagnetic phases (S-AF) as can be seen in the lower right part 
of the figure. A heuristic explanation of this phenomenon is given in the main text. 
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FIG. 16: Same phase diagram as in Fig.[^ this time obtained via Monte Carlo simulations on an 8^ -lattice. Like in the minimal model 
MF and Monte Carlo results agree quantitatively within an accuracy of 10% and less. The distinction between S and AC will become more 
apparent in FigtTTI 
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FIG. 17: Phase boundaries and orders of transitions as obtained via 
Monte Carlo simulation on a 8^ -lattice. We observe a mixture of 
both first and second order transitions depending on the particular 
values for the couplings. The symmetric phase is enclosed by ordered 
phases as already expected from the discussion of Section |lll| For 
further details on the simulation see the main text. 
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FIG. 18: Phase diagram of Fig.lTTlwith oriented trajectories (marked 
by arrows) used for the histograms of Figs.[l9lto|22l The curves are 
directed from the symmetric to the broken phases intersecting the 
critical lines vertically. We have chosen a representative subset from 
a total of twenty such curves which were analysed to determine the 
order of the transitions. 



tional freedom of choice between two (or more) states anti- 
aligned with a given one leads to an enormous degeneracy of 
the ground state in energy. As a consequence, entropy will be 
the sole judge deciding what is to be observed in a measure- 
ment. For the Z(3) Polyakov loop model both MF approxi- 
mation and the Monte Carlo simulations tell us that the most 
probable ground state corresponds to a preferred direction for 
the Polyakov loop on one of the two sub-lattices and an equal 
distribution of the two remaining directions on the other sub- 
lattice. Although we do not have an analytical justification for 
this statement the numerical evidence is compelling. Based 
on the latter, the sensitivity of ir to AF ordering can be ex- 
plained by a net expectation value which for the case of the 
3 -state Potts model is easily computed as 

£r = 0.5 zi + 0.25 Z2 + 0.25 ^3 7^ . (98) 

We conclude this discussion with an overview of the result- 
ing phases, their boundaries and the order of the transitions in 
between presented in Fig. [nl Our reasoning here is based 
on additional measurements carried out along parametrised 
curves A = X{s) as depicted in Fig. These simulations 
were exclusively focused on the order of the phase transition 
measuring the histograms of the observables L and M. As 
before all measurements were carried out on 8^ -lattices each 
trajectory X{s) being sampled with twenty points. For every 
such point Ag 10^ Monte Carlo sweeps were performed. To 
improve the statistics of the histograms we made use of all 
previously discussed symmetries (Z(3), complex conjugation 
and exchange of even and odd sub-lattice) by binning L or 
M together with their centre images zL, z'^L, zM, z'^M and 
their complex conjugates. This amounts to using each mea- 
sured value of L six times and of M even twelve times. In 



total we have recorded twenty such runs four of which are 
depicted in Fig.fTsl 

To illuminate the order of the transitions corresponding to 
the four directed line crossings displayed in Fig.[^we present 
six (out of twenty) histograms in Figs. [T9l to l22l which dis- 
play the distribution of the observable L respectively M in 
the complex plane. 

The first figure in this series. Fig. [T9| corresponds to the 
ferromagnetic transition, S-F. The histograms displayed may 
be viewed as a 'movie sequence' starting out in the symmetric 
phase with the distribution p{L) located at the origin. Re- 
calling the relation between probability distribution p{L) and 
the constraint effective potential lEHl . U{L) ex exp(— 
this situation corresponds to a unique minimum of U{L) at 
L = 0. As the couplings change three further maxima - 
which are Z(3) copies of each other - arise in addition to the 
one at the origin. Hence, we observe coexistence of ordered 
and disordered phases. The new maxima are separated from 
the original unique maximum at the origin by a finite amount. 
The associated discontinuity together with coexistence clearly 
shows that the transition is first order. 

In contrast to this the situation depicted in Fig.|20|is quite 
different. First of all the distribution is much more delocalised 
as compared to Fig. [191 implying ^ much flatter effective po- 
tential. Moreover, as soon as the maxima of the broken (or- 
dered) phase emerge the maximum at the origin has dissolved. 
Hence, there is no coexistence of different phases as in the pre- 
vious case but rather a continuous appearance of new maxima, 
branching away from the origin towards the corners of the fun- 
damental domain. These features characterise a second order 
phase transition. 

Analogous behaviour can be observed in Fig. I2TI with the 
maxima now moving in roughly the opposite directions in- 
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FIG. 19: Histogram of the observable L in the complex plane at a first order ferromagnetic phase transition. 
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FIG. 20: Histogram of the observable L in the complex plane at a second order ferromagnetic phase transition. 
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FIG. 21: Histogram of the observable L in the complex plane at a second order anti-centre phase transition. 
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FIG. 22: Histogram of the observable M in the complex plane at a second order anti-ferromagnetic phase transition. 
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dicating a transition to the AC phase. Still different is the 
S-AF transition of Fig. |22| By the same reasoning as given 
above this transition is of second order as well. However, the 
probability function emerging in the broken phase does not 
display the usual three-fold symmetry we have observed for 
the previous transitions. Although the figures seemingly look 
rotationally invariant one actually finds a six-fold degeneracy 
rather than a continuous one. The symmetry enhancement by 
a factor of two has a simple explanation. Note that Fig. |22| 
shows the probability distribution of the AF order parameter 
M which is defined in terms of even and odd sub-lattices. In- 
terchanging the latter accounts for the additional factor of two. 

Let us conclude this section with a short remark concerning 
Figs. |20| and Ell We have seen that our numerical data seem 
to indicate the possibility of second order transitions from the 
disordered phase S to F or AC ordering. This is at variance 
with the folklore claiming the absence of a Z(3) universal- 
ity class in d = 3, see e.g. It is not excluded that the 
observed second order characteristics near the tricritical point 
(S-F-AC) derive from artifacts due to our algorithms. At the 
moment, however, we cannot make a definite statement about 
this issue. Obviously, the puzzle deserves further investigation 
in the future. 



VI. SUMMARY AND OUTLOOK 

As outlined in the introduction the original motivation for 
Polyakov loop models lies in their status as effective theories 
for finite-temperature gluodynamics. However, our investiga- 
tion of the statistical mechanics involved should have made 
it clear that they are interesting in their own right. In sup- 
port of this statement we mention the 3 -fold Z(3) symmetry 
shared with the 3 -state Potts model, the nontrivial complex 
target space and the enormously rich phase structure resulting 
thereof. 

As we have seen, the latter can be understood qualitatively 
by simple energy-entropy considerations which predict a sym- 
metric phase S close to the origin in coupling space 'sur- 
rounded' by broken phases. These can exhibit ferromagnetic 
(F), anti-ferromagnetic (AF) or anti-centre (AC) ordering. 
This picture has been confirmed quantitatively both within a 
mean-field (MF) approximation and by extensive Monte Carlo 
simulations. Both methods find a first order transition for S-F 
and a second order one for S-AF. In order to capture the details 
of the phase diagram a total number of 8000 simulations (cor- 
responding to roughly 3000 hours of CPU time) was required. 
Naturally, this has added further refinements to our classical 
and MF analysis like precise values for critical couplings and 
exponents. 



The agreement between MF and Monte Carlo result is much 
better than naively expected for a statistical model in (i = 3. 
The S-F critical couplings agree to within 1%, the discontinu- 
ity M within 10% and the S-AF critical coupling within 20%. 
The precise agreement for the S-F transition suggests that the 
point where the S-F and S-AC critical lines merge is actually a 
tricritical point where two first order transitions merge into a 
second order one. It is known that the upper critical dimension 
for a tricritical point is three jiJIi^l so that in = 3 the MF 
approximation actually is exact apart from logarithmic correc- 
tions I44I1 . This observation is corroborated by the fact that the 
Polyakov loop model with two couplings is somewhat similar 
to the spin-one Blume-Capel ||45|,|4^ model which apart from 
an Ising term has an addition K 5?, Si = —1, 0, 1, with 
3 spin states. The model is known to have a tricritical point 
and is closely related to the 3-state Potts model ITtIi . The sec- 
ond order phase transition line is characterised by 3(i Ising 
critical exponents. Obviously, one should test the analogous 
exponents for our AC-F transition at the very left of our phase 
diagrams. 

One may also address the spatial localisation of the phases. 
On the 3(i lattice the coexisting phases should be separated by 
interfaces the tension of which can be measured (see e.g. Ilal 
and references therein). 

In the end one is of course interested in matching the ef- 
fective couplings to the underlying microscopic theory, i.e. 
finite-temperature Yang-Mills. The resulting curve Xij{(3) 
in the space of effective couplings should stay in the sym- 
metric and ferromagentic phases. We have already solved 
the matching-problem for the (simpler) case of SU{2) us- 
ing inverse Monte Carlo methods based on Schwinger-Dyson 
equations. For SU{3) this is considerably more difficult (at 
least technically) due to the increased complexity of the Haar 
measure. The problem is actually similar to one encountered 
in strong-coupling expansions where SU{2) group integrals 
can be performed analytically which is no longer true for 
SU (3). Nevertheless, we have been able to derive the rele- 
vant Schwinger-Dyson equations and hope to report on their 
applications in the near future. 
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